Main figures and tables
# Load original data extraction table
d = import((here("final_analyses", "output", "data", # file path
"respiratory-data.csv")), # file name
header = T)
d = d %>%
as_tibble() %>%
select(!starts_with("X"))
# Create different variables per subgroup
priors_simple_oxygen_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "simple oxygen",
outcome == "discharge"
)
priors_noninvasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "non-invasive ventilation",
outcome == "discharge"
)
priors_invasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "invasive ventilation",
outcome == "discharge"
)
recovery_simple_oxygen = d %>%
filter(trial == "RECOVERY",
outcome == "discharge",
oxygen == "simple oxygen")
recovery_noninvasive = d %>%
filter(trial == "RECOVERY",
outcome == "discharge",
oxygen == "non-invasive ventilation")
recovery_invasive = d %>%
filter(trial == "RECOVERY",
outcome == "discharge",
oxygen == "invasive ventilation")
# Load corticosteroids raw data
d_cortico = import(here("final_analyses", "data", # file path
"corticosteroids_subgroup_extracted-data.xlsx"))
d_cortico = clean_names(d_cortico)
priors_not_using =
d_cortico %>%
filter(
trial != "RECOVERY",
corticosteroids == 0,
outcome == "discharge"
)
priors_using =
d_cortico %>%
filter(
trial != "RECOVERY",
corticosteroids == 1,
outcome == "discharge"
)
recovery_not_using = d_cortico %>%
filter(trial == "RECOVERY",
outcome == "discharge",
corticosteroids == 0)
recovery_using = d_cortico %>%
filter(trial == "RECOVERY",
outcome == "discharge",
corticosteroids == 1)
# Set seed for reproducibility
set.seed = 123
Overview in priors and RECOVERY
respiratory_table =
d %>%
# only mortality
filter(outcome == "discharge") %>%
# change names to display them in the table
mutate(oxygen = case_when(
oxygen == "simple oxygen" ~ "Simple oxygen only",
oxygen == "non-invasive ventilation" ~ "Non-invasive ventilation",
oxygen == "invasive ventilation" ~ "Invasive mechanical ventilation",
)) %>%
# Remove any row with NA, eg, the ones with "Pooled" in oxygen
drop_na() %>%
# Create columns with Risk (%) for both control and tocilizumab
mutate(risk_control = 100*round(events_control/total_control, 2),
risk_toci = 100*round(events_toci/total_toci,2)) %>%
# Order columns
select(oxygen,
trial,
events_control, total_control, risk_control,
events_toci, total_toci, risk_toci) %>%
# Rename to match the next table
rename("Subgroup" = oxygen)
cortico_table =
d_cortico %>%
filter(outcome == "discharge") %>%
mutate(corticosteroids = case_when(
corticosteroids == 0 ~ "Not using corticosteroids",
corticosteroids == 1 ~ "Using corticosteroids"
)) %>%
mutate(risk_control = 100*round(events_control/total_control, 2),
risk_toci = 100*round(events_toci/total_toci,2)) %>%
select(corticosteroids,
trial,
events_control, total_control, risk_control,
events_toci, total_toci, risk_toci) %>%
rename("Subgroup" = corticosteroids)
# Put both tables togehter
both_tables = bind_rows(respiratory_table, cortico_table)
table1 = both_tables %>%
# Rename for final table
rename("Study" = trial,
"Events" = events_control,
"Total" = total_control,
"Risk (%)" = risk_control,
"Events " = events_toci,
"Total " = total_toci,
"Risk (%) " = risk_toci) %>%
# Define order for final table
mutate(Subgroup = factor(Subgroup,
levels = c("Not using corticosteroids",
"Using corticosteroids",
"Simple oxygen only",
"Non-invasive ventilation",
"Invasive mechanical ventilation")),
Study = factor(Study,
levels = c("RECOVERY",
"COVACTA",
"REMAP-CAP",
"CORIMUNO-19",
"Salvarini"))) %>%
# Re-arrange based on Study name
arrange(Study) %>%
# Gropup to create row tabs with gt()
group_by(Subgroup) %>%
# https://stackoverflow.com/questions/43832434/arrange-within-a-group-with-dplyr
arrange(Subgroup, .by_group = TRUE) %>%
### Transform into HTML table with gt()
gt() %>%
# Centralize columns
cols_align(
align = "center",
columns = everything()
) %>%
# Expand the size of first column
cols_width(
Study ~ px(250)
) %>%
# Add column tabs based on the column names defined above
tab_spanner(label = "Control",
columns = c("Events", "Total", "Risk (%)")) %>%
tab_spanner(label = "Tocilizumab",
columns = c("Events ", "Total ", "Risk (%) "))
table1
| Study | Control | Tocilizumab | ||||
|---|---|---|---|---|---|---|
| Events | Total | Risk (%) | Events | Total | Risk (%) | |
| Not using corticosteroids | ||||||
| RECOVERY | 168 | 367 | 46 | 162 | 357 | 45 |
| COVACTA | 35 | 65 | 54 | 124 | 188 | 66 |
| Using corticosteroids | ||||||
| RECOVERY | 873 | 1721 | 51 | 987 | 1664 | 59 |
| COVACTA | 36 | 79 | 46 | 42 | 106 | 40 |
| Simple oxygen only | ||||||
| RECOVERY | 635 | 933 | 68 | 697 | 935 | 75 |
| COVACTA | 35 | 44 | 80 | 66 | 78 | 85 |
| CORIMUNO-19 | 49 | 67 | 73 | 52 | 63 | 83 |
| Non-invasive ventilation | ||||||
| RECOVERY | 362 | 867 | 42 | 401 | 819 | 49 |
| COVACTA | 19 | 39 | 49 | 56 | 94 | 60 |
| Salvarini | 58 | 63 | 92 | 54 | 60 | 90 |
| Invasive mechanical ventilation | ||||||
| RECOVERY | 47 | 294 | 16 | 52 | 268 | 19 |
| COVACTA | 13 | 55 | 24 | 35 | 113 | 31 |
# table1 %>%
# gtsave(here("final_analyses", "output", "plots", "appendix", # File path
# "STable8_discharge_table1.png"))
## DEPRECATED
t = tibble(
"Subgroup" = c(
"Not using",
"Using",
"Simple oxygen only",
"Non-invasive ventilation",
"Invasive mechanical ventilation"
),
## Priors
"Number of studies" = c(
priors_not_using %>%
count() %>% pull(),
priors_using %>%
count() %>% pull(),
priors_simple_oxygen_raw %>%
count() %>% pull(),
priors_noninvasive_raw %>%
count() %>% pull(),
priors_invasive_raw %>%
count() %>% pull()
),
# Control arm
"Events " = c(
priors_not_using %>%
summarise(n = sum(events_control)) %>% pull(),
priors_using %>%
summarise(n = sum(events_control)) %>% pull(),
priors_simple_oxygen_raw %>%
summarise(n = sum(events_control)) %>% pull(),
priors_noninvasive_raw %>%
summarise(n = sum(events_control)) %>% pull(),
priors_invasive_raw %>%
summarise(n = sum(events_control)) %>% pull()
),
"Total " = c(
priors_not_using %>%
summarise(n = sum(total_control)) %>% pull(),
priors_using %>%
summarise(n = sum(total_control)) %>% pull(),
priors_simple_oxygen_raw %>%
summarise(n = sum(total_control)) %>% pull(),
priors_noninvasive_raw %>%
summarise(n = sum(total_control)) %>% pull(),
priors_invasive_raw %>%
summarise(n = sum(total_control)) %>% pull()
),
"Risk (%)" = c(
priors_not_using %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
priors_using %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
priors_simple_oxygen_raw %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
priors_noninvasive_raw %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
priors_invasive_raw %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull()
),
# Tocilizumab arm
"Events " = c(
priors_not_using %>%
summarise(n = sum(events_toci)) %>% pull(),
priors_using %>%
summarise(n = sum(events_toci)) %>% pull(),
priors_simple_oxygen_raw %>%
summarise(n = sum(events_toci)) %>% pull(),
priors_noninvasive_raw %>%
summarise(n = sum(events_toci)) %>% pull(),
priors_invasive_raw %>%
summarise(n = sum(events_toci)) %>% pull()
),
"Total " = c(
priors_not_using %>%
summarise(n = sum(total_toci)) %>% pull(),
priors_using %>%
summarise(n = sum(total_toci)) %>% pull(),
priors_simple_oxygen_raw %>%
summarise(n = sum(total_toci)) %>% pull(),
priors_noninvasive_raw %>%
summarise(n = sum(total_toci)) %>% pull(),
priors_invasive_raw %>%
summarise(n = sum(total_toci)) %>% pull()
),
"Risk (%) " = c(
priors_not_using %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
priors_using %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
priors_simple_oxygen_raw %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
priors_noninvasive_raw %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
priors_invasive_raw %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull()
),
## RECOVERY
# Control arm
"Events" = c(
recovery_not_using %>%
select(events_control) %>% pull(),
recovery_using %>%
select(events_control) %>% pull(),
recovery_simple_oxygen %>%
select(events_control) %>% pull(),
recovery_noninvasive %>%
select(events_control) %>% pull(),
recovery_invasive %>%
select(events_control) %>% pull()
),
"Total" = c(
recovery_not_using %>%
select(total_control) %>% pull(),
recovery_using %>%
select(total_control) %>% pull(),
recovery_simple_oxygen %>%
select(total_control) %>% pull(),
recovery_noninvasive %>%
select(total_control) %>% pull(),
recovery_invasive %>%
select(total_control) %>% pull()
),
"Risk (%) " = c(
recovery_not_using %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
recovery_using %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
recovery_simple_oxygen %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
recovery_noninvasive %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull(),
recovery_invasive %>%
summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
pull()
),
# Tocilizumab arm
"Events " = c(
recovery_not_using %>%
select(events_toci) %>% pull(),
recovery_using %>%
select(events_toci) %>% pull(),
recovery_simple_oxygen %>%
select(events_toci) %>% pull(),
recovery_noninvasive %>%
select(events_toci) %>% pull(),
recovery_invasive %>%
select(events_toci) %>% pull()
),
"Total " = c(
recovery_not_using %>%
select(total_toci) %>% pull(),
recovery_using %>%
select(total_toci) %>% pull(),
recovery_simple_oxygen %>%
select(total_toci) %>% pull(),
recovery_noninvasive %>%
select(total_toci) %>% pull(),
recovery_invasive %>%
select(total_toci) %>% pull()
),
"Risk (%) " = c(
recovery_not_using %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
recovery_using %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
recovery_simple_oxygen %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
recovery_noninvasive %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull(),
recovery_invasive %>%
summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
pull()
)
)
# Use kableExtra to have a nice table
# t1 = t %>%
# kbl(booktabs = T, align = 'c') %>%
# kable_styling(bootstrap_options = "striped", full_width = T) %>%
# column_spec(c(5,8,11,14), bold = T) %>%
# add_header_above(c(" " = 2,
# "Control" = 3,
# "Tocilizumab" = 3,
# "Control" = 3,
# "Tocilizumab" = 3)) %>%
# add_header_above(c(" " = 1,
# "Priors" = 7,
# "RECOVERY" = 6)) %>%
# pack_rows("Use of corticosteroids", 1, 2) %>%
# pack_rows("Respiratory support", 3, 5)
#
# t1
# t1 %>%
# save_kable(here("final_analyses", "output", "plots", "appendix", # File path
# "STable13_discharge_table1.pdf")) # File name
Prior and RECOVERY distributions
# Corticosteroids subgroups
# Create data frame
# Not using corticosteroids subgroup
distributions_not =
posteriors_not_using_discharge %>%
filter(type == "evidence-based") %>%
summarise(
RECOVERY_not = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior_not = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(RECOVERY_not:Prior_not)
distributions_not = bind_cols(
distributions_not,
# Bind corresponding posterior distribution data
samples_posteriors_not_using_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior_not" = `log(OR)`)
)
# Using corticosteroids subgroup
distributions_u =
posteriors_using_discharge %>%
filter(type == "evidence-based") %>%
summarise(
RECOVERY_u = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior_u = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(RECOVERY_u:Prior_u)
distributions_u = bind_cols(
distributions_u,
# Bind corresponding posterior distribution data
samples_posteriors_using_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior_u" = `log(OR)`)
)
# Bind altogether
distributions_cortico = bind_cols(
distributions_not,
distributions_u
)
arrow_labels = c("Tocilizumab Harm", "Tocilizumab Benefit")
xlim = c(0, 20)
xlab_df = data.frame(text = arrow_labels,
x = c(0.4,19),
y = c(0, 0),
hjust = c(0, 1))
a_small_amount = abs(xlim[1] - xlim[2])/20
null_line_at = 5
arrow_df = data.frame(id = c(1,2),
xstart = c(null_line_at - a_small_amount,
null_line_at + a_small_amount),
xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
y = c(1, 1))
arrows_plot = ggplot() +
geom_segment(data = arrow_df,
aes(x = .data$xstart,
xend = .data$xend,
y = .data$y,
yend = .data$y),
arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
geom_text(data = xlab_df,
aes(x = .data$x,
y = .data$y,
label = .data$text,
hjust = .data$hjust), size = 3.5) +
scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
scale_x_continuous(expand = c(0,0), limits = xlim) +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.box.background = element_rect(fill = "transparent"),
panel.border = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())
# Plot!
p1 = distributions_cortico %>%
# make it tidy for ggplot
pivot_longer(
cols = c(RECOVERY_not:Posterior_u),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Create new column to be able to use facet at the end
subgroup = case_when(
str_detect(dist, "_u") ~ "Using corticosteroids",
str_detect(dist, "_not") ~ "Not using corticosteroids"
),
# Set order of subgroups
subgroup = factor(subgroup,
levels = c(
"Not using corticosteroids",
"Using corticosteroids"
)
),
# Remove "_u" and "_not"
dist = case_when(
str_detect(dist, "RECOVERY") ~ "RECOVERY",
str_detect(dist, "Prior") ~ "Prior",
str_detect(dist, "Posterior") ~ "Posterior"
),
# Set order of distributions
dist = factor(dist,
levels = c(
"Posterior",
"RECOVERY",
"Prior"
)
)
) %>%
# exp() to convert log-odds ratio into odds ratio
ggplot(aes(
x = exp(samples), y = dist,
fill = dist
)) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(position = "dodge",
point_interval = median_qi, # Quantile interval
.width = c(0.95)) +
geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
scale_fill_manual(' ',
breaks=c("Prior","RECOVERY", "Posterior"),
values = c(
"gray70", # Prior
"#67B8D5", # RECOVERY
"#F8D08D" # Posterior
)) +
labs(
x = "\nOdds Ratio",
y = " "
) +
scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
coord_cartesian(xlim = c(0.25, 3),
clip = 'off') +
scale_y_discrete(expand = c(0, 0.5)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 11),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'right',
legend.text = element_text(size=12),
legend.key= element_blank(),
plot.margin = margin(20, 0, 70, 25)
) +
facet_wrap(~subgroup, ncol = 1)
## Respiratory support
# Create data frame
# Simple oxygen subgroup
distributions_s =
posteriors_simple_discharge %>%
filter(type == "evidence-based") %>%
summarise(
RECOVERY_s = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior_s = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(RECOVERY_s:Prior_s)
distributions_s = bind_cols(
distributions_s,
# Bind corresponding posterior distribution data
samples_posteriors_simple_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior_s" = `log(OR)`)
)
# Non-invasive ventilation
distributions_n =
posteriors_noninvasive_discharge %>%
filter(type == "evidence-based") %>%
summarise(
RECOVERY_n = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior_n = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(RECOVERY_n:Prior_n)
distributions_n = bind_cols(
distributions_n,
# Bind corresponding posterior distribution data
samples_posteriors_noninvasive_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior_n" = `log(OR)`)
)
# Invasive mechanical ventilation subgroup
distributions_i =
posteriors_invasive_discharge %>%
filter(type == "evidence-based") %>%
summarise(
RECOVERY_i = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior_i = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(RECOVERY_i:Prior_i)
distributions_i = bind_cols(
distributions_i,
# Bind corresponding posterior distribution data
samples_posteriors_invasive_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior_i" = `log(OR)`)
)
# Bind altogether
distributions_respiratory = bind_cols(
distributions_i,
distributions_n,
distributions_s
)
xlab_df = data.frame(text = arrow_labels,
x = c(0.4,19),
y = c(0, 0),
hjust = c(0, 1))
a_small_amount = abs(xlim[1] - xlim[2])/20
null_line_at = 5
arrow_df = data.frame(id = c(1,2),
xstart = c(null_line_at - a_small_amount,
null_line_at + a_small_amount),
xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
y = c(1, 1))
arrows_plot2 = ggplot() +
geom_segment(data = arrow_df,
aes(x = .data$xstart,
xend = .data$xend,
y = .data$y,
yend = .data$y),
arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
geom_text(data = xlab_df,
aes(x = .data$x,
y = .data$y,
label = .data$text,
hjust = .data$hjust), size = 3.5) +
scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
scale_x_continuous(expand = c(0,0), limits = xlim) +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.box.background = element_rect(fill = "transparent"),
panel.border = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())
# Plot!
p2 = distributions_respiratory %>%
# make it tidy for ggplot
pivot_longer(
cols = c(RECOVERY_i:Prior_s),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Create new column to be able to use facet at the end
subgroup = case_when(
str_detect(dist, "_s") ~ "Simple oxygen only",
str_detect(dist, "_n") ~ "Non-invasive ventilation",
str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
),
# Set order of subgroups
subgroup = factor(subgroup,
levels = c(
"Simple oxygen only",
"Non-invasive ventilation",
"Invasive mechanical ventilation"
)
),
# Remove "_i", "_n" and "_s"
dist = case_when(
str_detect(dist, "RECOVERY") ~ "RECOVERY",
str_detect(dist, "Prior") ~ "Prior",
str_detect(dist, "Posterior") ~ "Posterior"
),
# Set order of distributions
dist = factor(dist,
levels = c(
"Posterior",
"RECOVERY",
"Prior"
)
)
) %>%
# exp() to convert log-odds ratio into odds ratio
ggplot(aes(
x = exp(samples), y = dist,
fill = dist
)) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(position = "dodge",
point_interval = median_qi, # Quantile interval
.width = c(0.95)) +
scale_fill_manual(values = c(
"#F8D08D", # Posterior
"#67B8D5", # RECOVERY
"gray70" # Priors
)) +
geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
labs(
x = "\nOdds Ratio",
y = " "
) +
scale_x_continuous(breaks = seq(0.5, 3.5, 0.5)) +
coord_cartesian(xlim = c(0.8, 3.2)) +
scale_y_discrete(expand = c(0, 0.5)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 11),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 15, 70, 0)
) +
facet_wrap(~subgroup, ncol = 1)
p1_overlay = p1 + inset_element(arrows_plot,
ignore_tag = TRUE,
align_to = "full",
left = unit(1.5, 'cm'),
bottom = unit(0.5, 'cm'),
right = unit(11.3, 'cm'),
top = unit(2.5, 'cm'))
# Plot everything
p1_overlay +
p2 +
inset_element(arrows_plot2,
ignore_tag = TRUE,
align_to = "full",
left = unit(0, 'cm'),
bottom = unit(0.5, 'cm'),
right = unit(10.3, 'cm'),
top = unit(2.5, 'cm')) +
plot_annotation(tag_levels = 'A')
Point estimates depict the median and interval bars depict 95th quantiles.
# ggsave(width = 10,
# height = 6,
# here("final_analyses", "output", "plots", "appendix", # File path
# "SFigure13_discharge_figure1.pdf")) # File name
Here are the 95% quantile intervals in odds ratio for each distribution from the figure above.
t = distributions_cortico %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RECOVERY_not:Posterior_u),
names_to = "Subgroup", values_to = "samples"
) %>%
# Group by distribution so we can apply median_hdi
group_by(Subgroup) %>%
# Calculate the median and 95% HDI in log-odds ratio
ggdist::median_hdi(.width = 0.95) %>%
# Select only the relevant columns
select(Subgroup:.upper) %>%
# Exponentiate to transform log-odds ratio into odds ratio
mutate(across(samples:.upper, ~exp(.))) %>%
mutate(across(samples:.upper, ~round(., 2))) %>%
# Set order
mutate(Subgroup = factor(Subgroup,
levels = c(
"Prior_not", "RECOVERY_not", "Posterior_not",
"Prior_u", "RECOVERY_u", "Posterior_u"
))) %>%
arrange(Subgroup) %>%
rename("Median" = samples,
"Upper limit" = .upper,
"Lower limit" = .lower)
# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order
# create a vector with letters in the desired order
x = c("Prior_s", "RECOVERY_s", "Posterior_s",
"Prior_n", "RECOVERY_n", "Posterior_n",
"Prior_i", "RECOVERY_i", "Posterior_i")
t1 = distributions_respiratory %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RECOVERY_i:Posterior_s),
names_to = "Subgroup", values_to = "samples"
) %>%
# Group by distribution so we can apply median_hdi
group_by(Subgroup) %>%
# Calculate the median and 95% HDI in log-odds ratio
ggdist::median_hdi(.width = 0.95) %>%
# Select only the relevant columns
select(Subgroup:.upper) %>%
# Exponentiate to transform log-odds ratio into odds ratio
mutate(across(samples:.upper, ~exp(.))) %>%
mutate(across(samples:.upper, ~round(., 2))) %>%
# Define order
mutate(Subgroup = factor(Subgroup, levels = x)) %>%
arrange(Subgroup) %>%
rename("Median" = samples,
"Upper limit" = .upper,
"Lower limit" = .lower)
tfinal = bind_rows(t, t1) %>%
mutate(
Subgroup = case_when(
str_detect(Subgroup, "Prior") ~ "Prior",
str_detect(Subgroup, "RECOVERY") ~ "RECOVERY",
str_detect(Subgroup, "Posterior") ~ "Posterior"
)) %>%
rename("Distribution" = Subgroup) %>%
gt() %>%
cols_width(
Distribution ~ px(230)
) %>%
tab_row_group(
label = "Invasive mechanical ventilation",
rows = 13:15
) %>%
tab_row_group(
label = "Non-invasive ventilation",
rows = 10:12
) %>%
tab_row_group(
label = "Simple oxygen only",
rows = 7:9
) %>%
tab_row_group(
label = "Using corticosteroids",
rows = 4:6
) %>%
tab_row_group(
label = "Not using corticosteroids",
rows = 1:3
) %>%
cols_align(
align = "center",
columns = everything()
)
tfinal
| Distribution | Median | Lower limit | Upper limit |
|---|---|---|---|
| Not using corticosteroids | |||
| Prior | 1.66 | 0.94 | 2.96 |
| RECOVERY | 0.98 | 0.74 | 1.32 |
| Posterior | 1.10 | 0.84 | 1.42 |
| Using corticosteroids | |||
| Prior | 0.78 | 0.44 | 1.42 |
| RECOVERY | 1.42 | 1.24 | 1.62 |
| Posterior | 1.37 | 1.20 | 1.57 |
| Simple oxygen only | |||
| Prior | 1.59 | 0.83 | 2.94 |
| RECOVERY | 1.37 | 1.13 | 1.68 |
| Posterior | 1.39 | 1.15 | 1.68 |
| Non-invasive ventilation | |||
| Prior | 1.29 | 0.68 | 2.44 |
| RECOVERY | 1.34 | 1.10 | 1.62 |
| Posterior | 1.33 | 1.11 | 1.60 |
| Invasive mechanical ventilation | |||
| Prior | 1.45 | 0.70 | 3.03 |
| RECOVERY | 1.27 | 0.81 | 1.93 |
| Posterior | 1.31 | 0.90 | 1.92 |
# tfinal %>%
# gtsave(here("final_analyses", "output", "plots", "appendix", # File path
# "STable9_HDI_intervals_figureS13.png")) # File name
Posteriors distributions
# Create data frame
posterior_not_using =
samples_posteriors_not_using_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`RD`) %>%
rename("RD_not" = `RD`)
posterior_using =
samples_posteriors_using_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`RD`) %>%
rename("RD_u" = `RD`)
distributions_cortico_posterior =
bind_cols(posterior_not_using,
posterior_using)
# Prepare arrows under axis
# Adapted from https://github.com/rdboyes/forester/blob/master/R/forester.R
arrow_labels = c("Tocilizumab Benefit", "Tocilizumab Harm")
xlim = c(0, 20)
xlab_df = data.frame(text = arrow_labels,
x = c(2,19.5),
y = c(0, 0),
hjust = c(0, 1))
a_small_amount = abs(xlim[1] - xlim[2])/35
null_line_at = 14.5
arrow_df = data.frame(id = c(1,2),
xstart = c(null_line_at - a_small_amount,
null_line_at + a_small_amount),
xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
y = c(1, 1))
arrows_plot = ggplot() +
geom_segment(data = arrow_df,
aes(x = .data$xstart,
xend = .data$xend,
y = .data$y,
yend = .data$y),
arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
geom_text(data = xlab_df,
aes(x = .data$x,
y = .data$y,
label = .data$text,
hjust = .data$hjust), size = 3.5) +
scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
scale_x_continuous(expand = c(0,0), limits = xlim) +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.box.background = element_rect(fill = "transparent"),
panel.border = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())
p1 = distributions_cortico_posterior %>%
# make it tidy for ggplot
pivot_longer(
cols = c(RD_not:RD_u),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Create new column to be able to use facet at the end
subgroup = case_when(
str_detect(dist, "_u") ~ "Using corticosteroids",
str_detect(dist, "_not") ~ "Not using corticosteroids"
),
# Set order of subgroups
subgroup = factor(subgroup,
levels = c(
"Not using corticosteroids",
"Using corticosteroids"
)
)) %>%
# https://mjskay.github.io/ggdist/articles/slabinterval.html
ggplot(aes(
# Multiply by 100 to report in percentage
x = 100 * samples, y = dist,
fill = subgroup,
# https://github.com/mjskay/ggdist/issues/71
fill_ramp = stat(x < 0))
) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(point_interval = median_hdi, # Highest density interval
.width = c(0.8, 0.95),
# https://github.com/mjskay/ggdist/issues/70
interval_size_range = c(1.1,1.9)) +
scale_fill_ramp_discrete(from = "gray85", range = c(0,1)) +
scale_fill_manual(values = c("#9EA45B", "#A65041")) +
geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
labs(
x = "\nRisk Difference (%)",
y = " "
) +
scale_x_continuous(breaks = seq(from = -15, to = 5, 2.5)) +
coord_cartesian(xlim = c(-15, 5)) +
scale_y_discrete(expand = c(0, 0.5)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 11),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(10, 10, 50, 10)
) +
facet_wrap(~subgroup, ncol = 1, scales = "free_y")
# Create data frame
# Simple oxygen only
posterior_s =
samples_posteriors_simple_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`RD`) %>%
rename("RD_s" = `RD`)
# Non-invasive ventilation
posterior_n =
samples_posteriors_noninvasive_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`RD`) %>%
rename("RD_n" = `RD`)
# Invasive mechanical ventilation
posterior_i =
samples_posteriors_invasive_discharge %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`RD`) %>%
rename("RD_i" = `RD`)
# Bind altogether
distributions_respiratory_posterior = bind_cols(
posterior_i,
posterior_n,
posterior_s
)
p3 = distributions_respiratory_posterior %>%
# make it tidy for ggplot
pivot_longer(
cols = c(RD_i:RD_s),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Create new column to be able to use facet at the end
subgroup = case_when(
str_detect(dist, "_s") ~ "Simple oxygen only",
str_detect(dist, "_n") ~ "Non-invasive ventilation",
str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
),
# Set order of subgroups
subgroup = factor(subgroup,
levels = c(
"Simple oxygen only",
"Non-invasive ventilation",
"Invasive mechanical ventilation"
))
) %>%
# Plot!
ggplot(aes(
# Multiply by 100 to report in percentage
x = 100 * samples, y = dist,
fill = subgroup,
# https://github.com/mjskay/ggdist/issues/71
fill_ramp = stat(x < 0))
) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(point_interval = median_hdi, # Highest density interval
.width = c(0.8, 0.95),
# https://github.com/mjskay/ggdist/issues/70
interval_size_range = c(1.1,1.9)) +
scale_fill_ramp_discrete(from = "gray85", range = c(0,1)) +
scale_fill_manual(values = c("#9D75B1", # Simple oxygen
"#CE9F5B", # Non-invasive
"#5891BF") # Invasive
) +
geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
labs(
x = "\nRisk Difference (%)",
y = " "
) +
scale_x_continuous(breaks = seq(from = -15, to = 5, 2.5)) +
coord_cartesian(xlim = c(-15, 5)) +
scale_y_discrete(expand = c(0, 0.5)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 11),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 10, 50, 10)
) +
facet_wrap(~ subgroup, ncol = 1, scales = "free_y")
p4 = distributions_respiratory_posterior %>%
# make it tidy for ggplot
pivot_longer(
cols = c(RD_i:RD_s),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Create new column to be able to use facet at the end
subgroup = case_when(
str_detect(dist, "_s") ~ "Simple oxygen only",
str_detect(dist, "_n") ~ "Non-invasive ventilation",
str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
),
# Set order of subgroups
subgroup = factor(subgroup,
levels = c(
"Simple oxygen only",
"Non-invasive ventilation",
"Invasive mechanical ventilation"
))
) %>%
# Multiply by 100 to report in percentage
ggplot(aes(100*samples, colour = subgroup)) +
stat_ecdf(geom = "step", size = 1.2) +
scale_color_manual(' ',
values = c("#9D75B1", # Simple oxygen
"#CE9F5B", # Non-invasive
"#5891BF") # Invasive
) +
geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
labs(
x = "\nRisk Difference (%)",
y = "Probability RD \u2264 X\n"
) +
scale_x_continuous(
breaks = seq(from = -15, to = 10, 1),
limits = c(-15, 10)) +
scale_y_continuous(
breaks = seq(from = 0, to = 1, .1),
expand = c(0, .03)
) +
theme(
axis.ticks.x = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
panel.grid.major.y = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(10, 10, 50, 10)
) +
coord_cartesian(x = c(-12.5,0.1))
p1_overlay = p1 + inset_element(arrows_plot,
ignore_tag = TRUE,
align_to = "full",
left = unit(0.65, 'cm'),
bottom = unit(0, 'cm'),
right = unit(12.3, 'cm'),
top = unit(1.9, 'cm'))
# Plot everything
(p1_overlay + p2) / (p3 +
inset_element(arrows_plot,
ignore_tag = TRUE,
align_to = "full",
left = unit(0.65, 'cm'),
bottom = unit(0, 'cm'),
right = unit(12.3, 'cm'),
top = unit(1.9, 'cm')) +
p4) +
plot_annotation(tag_levels = 'A')
Panels A and C: Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals. Panels B and D: Cumulative posterior distributions correspond to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
# ggsave(width = 10,
# height = 8.5,
# #https://stackoverflow.com/questions/12768176/unicode-characters-in-ggplot2-pdf-output
# device=cairo_pdf,
# here("final_analyses", "output", "plots", "appendix", # File path
# "SFigure14_discharge_figure2.pdf")) # File name
multiply100 = function(x){x*100}
t1 = distributions_cortico_posterior %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RD_not:RD_u),
names_to = "Subgroup", values_to = "samples"
) %>%
group_by(Subgroup) %>%
summarise("Pr(\u2264 0%)" = (sum( samples <= 0))/ n(),
"Pr(\u2264 -1%)" = (sum( samples <= -0.01))/ n(),
"Pr(\u2264 -2%)" = (sum( samples <= -0.02))/ n(),
"Pr(\u2264 -5%)" = (sum( samples <= -0.05))/ n()) %>%
mutate(across(-Subgroup, ~multiply100(.))) %>%
mutate(across(-Subgroup, ~round(.,1))) %>%
mutate(
Subgroup = case_when(
str_detect(Subgroup, "_u") ~ "Using",
str_detect(Subgroup, "_not") ~ "Not using"
))
multiply100 = function(x){x*100}
t2 = distributions_cortico_posterior %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RD_not:RD_u),
names_to = "Subgroup", values_to = "samples"
) %>%
# Group by distribution so we can apply median_hdi
group_by(Subgroup) %>%
# Calculate the median and 95% quantile interval in log-odds ratio
ggdist::median_hdi(.width = 0.95) %>%
# Select only the relevant columns
select(Subgroup:.upper) %>%
mutate(across(samples:.upper, ~multiply100(.))) %>%
mutate(across(samples:.upper, ~round(., 2))) %>%
rename("Median" = samples,
"Upper limit" = .upper,
"Lower limit" = .lower) %>%
mutate(
Subgroup = case_when(
str_detect(Subgroup, "_u") ~ "Using",
str_detect(Subgroup, "_not") ~ "Not using"
))
multiply100 = function(x){x*100}
t3 = distributions_respiratory_posterior %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RD_i:RD_s),
names_to = "Subgroup", values_to = "samples"
) %>%
group_by(Subgroup) %>%
summarise("Pr(\u2264 0%)" = (sum( samples <= 0))/ n(),
"Pr(\u2264 -1%)" = (sum( samples <= -0.01))/ n(),
"Pr(\u2264 -2%)" = (sum( samples <= -0.02))/ n(),
"Pr(\u2264 -5%)" = (sum( samples <= -0.05))/ n()) %>%
mutate(across(-Subgroup, ~multiply100(.))) %>%
mutate(across(-Subgroup, ~round(.,1))) %>%
arrange(desc(Subgroup)) %>%
mutate(
Subgroup = case_when(
str_detect(Subgroup, "_s") ~ "Simple oxygen only",
str_detect(Subgroup, "_n") ~ "Non-invasive ventilation",
str_detect(Subgroup, "_i") ~ "Invasive mechanical ventilation"
))
multiply100 = function(x){x*100}
t4 = distributions_respiratory_posterior %>%
# make it tidy to be able to use group_by()
pivot_longer(
cols = c(RD_i:RD_s),
names_to = "Subgroup", values_to = "samples"
) %>%
# Group by distribution so we can apply median_hdi
group_by(Subgroup) %>%
# Calculate the median and 95% quantile interval in log-odds ratio
ggdist::median_hdi(.width = 0.95) %>%
# Select only the relevant columns
select(Subgroup:.upper) %>%
mutate(across(samples:.upper, ~multiply100(.))) %>%
mutate(across(samples:.upper, ~round(., 2))) %>%
rename("Median" = samples,
"Upper limit" = .upper,
"Lower limit" = .lower) %>%
arrange(desc(Subgroup)) %>%
mutate(
Subgroup = case_when(
str_detect(Subgroup, "_s") ~ "Simple oxygen only",
str_detect(Subgroup, "_n") ~ "Non-invasive ventilation",
str_detect(Subgroup, "_i") ~ "Invasive mechanical ventilation"
))
Here are the posterior probabilities for each posterior distribution.
bind_rows(t1,t3) %>%
kbl(booktabs = T, align = 'c') %>%
kable_styling(bootstrap_options = "striped", full_width = F) %>%
pack_rows("Use of corticosteroids", 1, 2) %>%
pack_rows("Respiratory support", 3, 5) %>%
add_header_above(c(" " = 1,
"Probability of Risk Difference \u2265 X%" = 4))
| Subgroup | Pr(≤ 0%) | Pr(≤ -1%) | Pr(≤ -2%) | Pr(≤ -5%) |
|---|---|---|---|---|
| Use of corticosteroids | ||||
| Not using | 75.6 | 65.2 | 53.6 | 20.8 |
| Using | 100.0 | 100.0 | 100.0 | 95.9 |
| Respiratory support | ||||
| Simple oxygen only | 100.0 | 99.8 | 99.3 | 82.2 |
| Non-invasive ventilation | 99.9 | 99.6 | 98.6 | 81.8 |
| Invasive mechanical ventilation | 92.1 | 84.8 | 74.8 | 36.8 |
Here are the 95% HDI in risk difference of each posterior distribution.
bind_rows(t2,t4) %>%
kbl(booktabs = T, align = 'c') %>%
kable_styling(bootstrap_options = "striped", full_width = F) %>%
pack_rows("Use of corticosteroids", 1, 2) %>%
pack_rows("Respiratory support", 3, 5) %>%
add_header_above(c(" " = 1,
"95% Highest Density Interval" = 3))
| Subgroup | Median | Lower limit | Upper limit |
|---|---|---|---|
| Use of corticosteroids | |||
| Not using | -2.30 | -8.77 | 4.16 |
| Using | -7.87 | -11.04 | -4.63 |
| Respiratory support | |||
| Simple oxygen only | -6.75 | -10.27 | -3.08 |
| Non-invasive ventilation | -7.13 | -11.74 | -2.56 |
| Invasive mechanical ventilation | -3.96 | -10.33 | 1.58 |